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Abstract 

Orbital and rotational variations perturb the latitudinal and seasonal pattern of incident 
solar radiation, producing major climatic change on time scales of 10 4 -10 6 years. The orbital 
variations are oblivious to internal structure and processes, but the rotational variations are not. 
The intent of this article is to describe a program of investigation whose objective would be to 
explore and quantify three aspects of orbital, rotational and climatic interactions. An important 
premise of this investigation is the synergism between geodynamics and paleoclimate. Better 
geophysical models of precessional dynamics are needed in order to accurately reconstruct the 
radiative input to climate models. Some of the paleoclimate proxy records contain information 
relevant to solid Earth processes, on time scales which are difficult to constrain otherwise. 

Specific mechanisms which will be addressed include: 

• climatic consequences of deglacial polar motion, and 

• precessional and climatic consequences of 

— glacially induced perturbations in the gravitational oblateness, and 

— partial decoupling of the mantle and core. 

The approach entails constructing theoretical models of the rotational, deformational, ra- 
diative and climatic response of the Earth to known orbital perturbations, and comparing these 
with extensive records of paleoclimate proxy data. Several of the mechanisms of interest may 
participate in previously unrecognized feed-back loops in the climate dynamics system. A new 
algorithm for estimating climatically diagnostic locations and seasons from the paleoclimate 
time series is proposed. 
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INTRODUCTION 


It is widely recognized that mass redistributions associated with climatic change (glaciations) 
are an important source of crustal deformation and geodynamic change. It is much less widely 
appreciated that rates, phases and amplitudes of deformation of the deep interior of the Earth 
can influence climate. The objective of this investigation is to better characterize four aspects of 
this geodynamic contribution to global climatic change. The common theme among them has two 
threads: internal mass redistributions influence the rotational dynamics of the Earth, and changes 
in orbital and rotational parameters influence the latitudinal and seasonal pattern of insolation. 
Previous attempts to account for astronomically forced climatic change have usually only considered 
extremely simplisitic models for the response of the Earth to external torques and surface loads. 

The focus of the proposed investigation will consist of three specific aspects of interaction between 
internal mass flow and rotational dynamics, and a new algorithm for comparing paleoclimatic proxy 
data with astronomically forced variations in insolation patterns. 

The latitudinal and seasonal pattern of incident solar radiation depends on the eccentricity of the 
Earth’s orbit and the orientation of the spin axis relative to both the orbit normal and the aspidal 
line. Unit vectors s and n characterize the directions of the spin axis and orbit normal, respectively. 
Two angles completely characterize the relative orientation of the spin axis. The obliquity c is simply 
the angle between the orbit normal and the spin axis 

e = cos'^n • s) (1) 

The ascending node of the orbit plane on the instantaneous equator plane has an orientation given 
by (s x n), and the longitude of perihelion zv is just the angle in the orbit plane from that node to 
perihelion. It is widely appreciated that secular variations in these three parameters (e,e,r?) produce 
major climatic change (Hays et al., 1976; Berger et ah, 1984). In fact, spectral analyses of long, high 
resolution marine sediment isotopic records show significant variance at periods near 100 kyr, 41 kyr 
and 19-23 kyr, which are generally attributed to spectral lines in the radiative forcing fluctuations 
associated with e, e and e sin(iu), respectively. 

The causes and effects of the orbital changes are quite well understood. Gravitational interactions 
with the other planets cause the shape and orientation of the orbit to change on time scales of 
10 4 -10 6 years. The inclination I and nodal longitude Q determine the orientation of the orbit plane. 
The eccentricity e and perihelic longitude oj determine the shape of the orbit and its orientation 
within the plane. Note that w is measured from an inertially fixed direction, rather than the moving 
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node as is the case for w. The secular evolution of the orbital element pairs (I,fi) and (e,w) can be 
conveniently represented in terms of Poisson series 


p zz sin(I) sin(fi) 

= £ Nj sin(sj t + gj) 

(2) 

q = sirt(I) cos(Q) 

= £ Nj cos(sj t + gj) 

h — e sin(u/) = 

£ Mj sin(rj t + fj) 

(3) 

k = e cos(w) = 

£ Mj cos(rj t fj) 


In the lowest order solution, there are as many frequencies rj and Sjas there are planets. However, 
the frequencies rj and s ; are characteristic modal frequencies (eigenvalues) of the coupled system 
of oscillators and are not each uniquely associated with a particular planet (Milani, 1988). The 
frequencies rj are all positive, indicating that the perihelia advance. In the lowest order solution, 
the apsidal rates are all in the interval (0.667 < rj < 28.221 arcsec/year). The corresponding 
periods are 45.92 kyr to 1.943 Myr. One of the frequencies Sj is zero, and all the others are negative, 
indicating that the nodes regress. In the lowest order solution, the non-zero nodal rates are all in 
the interval (0.692 < Sj < 26.330 arcsec/year). The corresponding periods are 49.22 kyr to 1.873 
Myr. In higher order solutions, variations in (e,w) become coupled to variations in (1,0), but the 
solutions can still be cast in terms of Poisson series like Equations (2) and (3). 

Laskar (1988) has recently published a secular variation theory which is complete to fifth order in 
eccentricity and inclination. Agreement between this secular variation model and strictly numerical 
computations (Richardson and Walker, 1989; Quinn et al., 1991) is much better than for any previous 
analytical model. The inclination and eccentricity series for Earth each contain 80 distinct terms. 
Figures (1) and (2) illustrate the spectra and corresponding histories of variation in orbital inclination 
and eccentricity for the past 2 10 6 years. 

In computing these secular orbital variations, the Earth, Moon and planets can all be treated as 
point masses. No internal structure or processes are relevant to orbital evolution. The physics of 
the process is simple, and well understood, though development of proper mathematical tools to 
represent the long term evolution remains an area of active research (Laskar, 1988, 1990; Quinn et 
al., 1991). On the other hand, the rotational evolution does depend rather sensitively on various 
aspects of the structure and dynamics of the interior. 

Lunar and solar gravitational torques acting on the oblate figure of the Earth cause the spin axis 
s to precess about the instantaneous orbit normal n. If the Earth is considered to be a rigid body, 
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the evolution of the spin axis orientation is given by 


ds/dt = a(n • s)(s x n) 


( 4 ) 


where 


%(C - A) Gm t (1 „ 2 




2Cn ^ K1,n ^ 

is a scalar rate factor which depends on intrinsic properties of the Earth, such as polar and 
equatorial moments of inertia (G,A) and rotation rate n, and on extrinsic influences, such as masses 
m, orbital inclinations I, and semiminor axes b, of the Moon and Sun. The solar and lunar torques 


together produce a precession of the spin axis of the Earth at a rate of c*(n ■ s) = 50.38 arcsec/year 
(Kinoshita, 1977, Williams et al., 1991). 

Once the present spin axis direction s is known and orbital element histories are given via Equations 
(2) and (3), an obliquity history can be constructed from equation (4) in two different ways. The 
linear perturbation approach (Miskovic, 1931, Sharaf and Boudnikova, 1967; Vernekar, 1972, 1977; 
Ward, 1974; Berger, 1976) involves deriving coefficients of a trigonometric series, similar to Equations 
(2) and (3), which yield the obliquity and longitude of perihelion directly as functions of time. An 
alternative is to apply standard numerical algorithms for solving initial value problems to generate 
a vector time series s(t) and then compute the obliquity and longitude of perihelion directly (Ward, 
1979, Laskar, 1986; Bills, 1990b). Figure (3a) shows the spectrum of obliquity variations, which 
in the linear perturbation model is simply obtained from the inclination spectrum by shifting each 
frequency s j by the luni-solar precession rate (s = 50.38 arcsec/year) and multiplying each amplitude 
Nj by the spectral admittance 


Fj — s j/( s j T s ) (6) 

Figure (3b) shows the numerically integrated obliquity history of Laskar (1986,1988) for the last 
2 Ma, using the rotational theory of Kinoshita (1977). The difference between the numerical and 
linear perturbation solution never exceeds 0.06 degree over that interval. It is clear that the linear 
perturbation solution gives a very adequate representation of the spin precession. 

The influence of orbital and rotational variations on climate is operative through perturbations 
in the latitudinal and seasonal pattern of insolation. The diurnal average intensity of radiation at a 
point is inversely proportional to the squared solar distance and directly proportional to the diurnal 
average rectified solar direction cosine 


F — ( a / r ) 2 (ll u • u 5 ||) 


(7) 
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where a and r are mean and instantaneous solar distance and u and u 5 are unit vectors from 
the center of the Earth to the surface point of interest and the sub-solar point, respectively. The 
insolation pattern, a s a function of latitude 6 and mean anomaly M, can be readily computed 
once values are specified for the orbital and rotational parameters e, e, and w (Hargreaves, 1895; 
Milankovitch, 1920; Vernekar, 1972, 1977; Ward, 1974). Figure (4) illustrates that pattern for 
the present orbital and rotational configuration. This pattern can also be written in terms of a 
Fourier-Legendre series (Hargreaves, 1895; North and Coakley, 1979; Taylor, 1984; Bills, 1992) 

F(fi, M-,e,e,ix) = P n (fi)J2exp(ip M) F„ lP (e,e,w) (8) 

where ji = cos(0), and P n is a Legendre polynomial. The number of terms in the Fourier summa- 
tion required to obtain a good representation of the seasonal pattern is greater in the polar regions 
than in the tropics and mid-latitudes. The primary difficulty in the polar regions is reproducing 
the abrupt change in slope of the insolation curve at times of transition to continual darkness or 
continual light. It is also true that the polar regions place the greatest demands on the Legendre 
summation, since the spatial pattern also has a discontinuous first derivative at the latitude where 
the transition occurs to continual darkness or light. 

A significant fraction of the recent work on comparing paleoclimate proxy records to astronomically 
forced insolation changes has been based on the published insolation curves of Berger (1978, 1991). 
Common practice is to compare computed variations in some particular aspect of the seasonal and 
latitudinal insolation pattern (July insolation at 65° N is a particularly frequent choice) with an 
observational record of some climatic indicator (<5 18 0 variations versus age (depth) in a marine 
sediment core, for example). Comparisons of this sort enable estimates of amplitude, phase and 
coherence of climatic response to radiative forcing, and have unequivocally demonstrated that orbital 
and rotational variations are a dominant cause of climatic change on time scales of 10 4 -10 6 years. 

Despite obvious successes, several problems remain in the general methodology. For example, the 
orbit and rotation are both assumed to be perfect clocks and, as a result, data time series are often 
” tuned” to the astronomical time scale (Martinson et ah, 1982; Shackleton et al., 1990; Hilgren, 
1991). However, as we shall see in the next two sections, the rotational variations in particular 
do not keep very good time. Also, a number of significant observations remain without adequate 
explanation. One of the most perplexing of these is the mid- Pleistocene climatic switch in dominant 
oscillation frequency. The 41 kyr oscillation was dominant throughout the early Pleistocene, and 
the 100 kyr oscillation has been dominant since the Brunhes-Matuyama magnetic polarity reversal 
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(or thereabouts), whereas the standard July 65° N insolation curve shows no discernible change in 
spectral composition over the entire interval. 

The intent of this article is to point out the need for a critical examination of this standard scenario, 
with the objective of improving the geodynamic component of the model enough to provide better 
radiative forcing time series to the paleoclimate community, use paleoclimate proxy data records to 
calibrate solid earth responses to applied loads and torques, and explore potential feed-back loops . 

Two of the topics of study are perturbations to the simple precession model presented in Equation 
(4). The first topic concerns time variations in the precession rate occasioned by glacial mass 
transport and resulting changes in the difference between polar and equatorial moments of inertia 
(C-A). These fractional variations can likely reach 1% in magnitude and are fully competitive with 
changes in orbital eccentricity in terms of their effect on instantaneous precession rate. The second 
topic is the effect of differential precession in the deep interior, which is governed by the (poorly 
known) characteristics of the inertial and dissipative coupling torques which attempt to keep the 
rotation axes of the mantle and core aligned. 

The third topic is the potentially significant geodynamic feed-back loop associated with deglacial 
mass transport, polar motion and ensuing perturbations to radiative equilibrium temperature pat- 
terns. The asymmetrical disposition of major ice sheets and ocean basins relative to the rotation axis 
implies that during growth or decay of these ice sheets the geographic location of the principal axis 
of greatest inertia will shift, ff the hydrospheric and cryospheric changes are appreciably more rapid 
than can be compensated by asthenospheric flow, the rotation axis will shift, possibly by as much 
as 1°. The primary climatic consequence is a shifting of the geographic distribution of continental 
and oceanic regions (which differ by a factor of 60 in heat capacity). The effect of a 1° change in 
the geographic position of the pole is different (but not necessarily less significant) than a 1° change 
in obliquity. 

The fourth and final topic is a new modelling strategy in which climate proxy data records are 
used to estimate linear combinations of insolation pattern Fourier- Legendre coefficients which best 
duplicate the observed variations. This approach allows the data variance to be partitioned into 
global versus regional effects, and can distinguish between responses to annual average insolation 
and those due to seasonal cycle fluctuations. 
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PRECESSION AL DYNAMICS WITH VARIABLE RATE FACTOR 


Statement of Problem 

All but the most recent reconstructions of the radiative forcing input to paleoclimate models have 
assumed that both the orbital and rotational dynamics could be readily and accurately reconstructed 
from their present configurations, via the simple analyses mentioned in the introduction. These 
expectations seem well founded in the case of orbital evolution, though the possibility of chaotic 
dynamics in the inner solar system (Laskar, 1990; Laskar et ah, 1992) does seem to preclude confident 
extrapolation beyond 10 7 years. However, there are a number of processes, working in different 
locations and at different rates, which all serve to compound the difficulty of accurately computing 
the spin precessionai evolution. 

An important aspect of this problem is the synergism between geodynamics and paleoclimate. 
Better geophysical models are needed in order to accurately reconstruct the radiative input to 
climate models. Some of the paleoclimate proxy records contain information relevant to solid Earth 
processes on time scales which are difficult to constrain otherwise. 

On the longest time scales of interest (10 7 -10 9 years) the limiting uncertainty is variability in the 
tidal transfer of angular momentum from the rotation of the Earth to the orbit of the Moon. At 
present, these tidal torques are increasing the length of the day by 22.5 10” 6 sec/year and increasing 
the size of the lunar orbit by 3.88 cm/year (Cazenave and Daillet, 1981; Christodoulidis et ah, 
1988). Berger et ah (1989) have made a useful first step towards including this effect in climatic 
time series. They computed the change in the major precession and obliquity frequencies due to 
lunar tidal evolution, assuming that the present rate of tidal energy dissipation is representative of 
the past 500 Myr. However, the present rates are considerably higher than the long term average 
(Hansen, 1982), largely due to a near resonance between sloshing modes of ocean basins and the 
diurnal and semidiurnal tidal periods (Platzman et ah, 1981), and apparently compounded by a 
contribution from shallow seas (Wunsch, 1986; Dickman and Preisig, 1986). Sedimentary records 
which constrain lunar orbital evolution show some promise of resolving this problem (Olsen, 1986; 
Williams, 1989a, b; Herbert and D’Hondt, 1990), but the situation is definitely more complex than 
is suggested by Berger et ah, (1989). 

Another parameter which can vary, on rather shorter time scales and in an equally irregular 
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fashion, is the gravitational oblateness of the Earth (C-A)/C. Thomson (1990) has recently made 
three important contributions to the understanding of this source of variability, first, he pointed 
out that, mass redistribution associated with major glaciations and compensating subsidence and 
crustal deformations (Le Treut and Ghil, 1983; Wu and Peltier, 1984) can cause fractional changes 
in oblateness of order 10 -3 -10~ 2 . Second, he showed that high resolution spectral analyses of 
several climatic time series appear to indicate fluctuations of the luni-solar precession rate of this 
magnitude, and with a dominant period near 100 kyr. Finally, Thomson pointed out that the 
best fit. to the paleoclimate proxy data was obtained using a mean lunisolar precession rate 0.6 
arcsec/year less than the present observed value. He notes that the resulting value would correspond 
rather closely with that expected for a hydrostatic flattening (Nakiboglu, 1982). If these important 
results are corroborated, they will demonstrate that important feed-back loops exist in the orbital- 
rotationai-climatic interactions system, further ” up-stream” in the presumed causal chain than has 
been previously recognized. 

Approach 

The research design for this segment of the proposed investigation will address several issues. The 
primary focus will be an attempt to resolve two related questions. What are the precessional and 
paleoclimatic consequences of small (0.00 1-0.0 1 ) variations in oblateness (C-A)/C, over time scales 
of 1 0 3 - 1 0 6 years? Can such variations can be confidently inferred from paleoclimatic proxy data? 
These questions represent a forward modeling problem and a coupled inverse problem, respectively. 

The first question is the easier to answer. Equation (4) describes the variations in orientation of 
the spin axis, as viewed in an inertial reference frame. However, since we are primarily interested in 
the orientation of the spin axis relative to orbit normal (obliquity) and the apsidal line (longitude 
of perihelion), the analysis will be made easier if we first transform to a coordinate system fixed in 
the orbit plane. If the rotation matrix is denoted by A, the transformed equation takes the form 
(Ward, 1974; Bills, 1990a) 

ds/dt = {d.A/dt ,4 -1 }s -f a(n ■ s)(s x n) (9) 

where 

{dA/dt A~ 1 } - B dl/dt + C dW/dt (10) 
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0 0 0 


( 11 ) 


B= 0 0 1 

0-10 

0 cos(7) — sin (7) 

C = — cos(7) 0 0 (12) 

sin(7) 0 0 

Now define two complex quantities: P = sin(I) e in , which represents components of the orbit 
normal on the invariable plane, and S = sin(e) e*’ 1 ', which represents components of the spin vector 
on the orbital plane. In this new notation, Equation (9) can be rewritten in the form 

dS/dt + ia S = ib dPjdt (13) 

where 

a = acos(e) + cos(I) dQ/dt (14) 

b — cos(c)e~ lil (15) 

The complete solution to a nonhomogeneous linear differential equation consists of both ”free” 
and ” forced 1 ’ modes of oscillation. The free modes, in this case, correspond to spin precession with 
the orbit plane fixed, and the forced modes correspond to motions of the spin axis as it attempts to 
precess about the instantaneous orbit normal, while the orbit normal itself is precessing. The forced 
modes make first order contributions to both e and zu, whereas the free modes are only second order 
for e but are first order for w. As a result, in the standard linear series solutions (Berger, 1978; 
Ward, 1974), the obliquity terms include only the forced response to changes in (I,fi), whereas the 
nodal longitude terms include both forced modes from (I,fi), and free modes with variable (e,u>). 

The result of variable oblateness will have exactly the same qualitative effect on spin precession as 
does a change in orbital eccentricity. Both effect the free modes only. It will thus be rather difficult 
to confidently distinguish oblateness variations from unmodeled eccentricity variations. However, 
that does not effect the forward modeling aspect of the problem. To the extent that oblateness 
variations occur, their effect should be included in the astronomical forcing to climate models. 

An iterative approach to the problem seems promising. Orbital variations are uneffected by 
oblateness and need not directly concern us. As a first step, standard rotational variations (including 
eccentricity, but neglecting oblateness variations) will be used to generate radiative input to an 
energy balance climate model (North et al., 1981; Short et al., 1991), with a coupled ice-sheet model 
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(DeBlonde and Peltier, 1990). Time variations in oblateness can be simply estimated from the 
surface load and internal compensation. The resulting oblateness history is then fed back into the 
rotational calculation and the entire process is repeated. The inner-most loops of this algorithm are 
somewhat similar to the model of Peltier (1982). However, he did not allow the mass loading to 
influence the radiative forcing. 

Even fairly modest changes in oblateness are rotationally significant. For comparison, the eccen- 
tricity perturbation influences precession rate via the factor (a/b) 3 , where a and b are semimajor 
and semiminor axes, respectively. This amounts to (l-e 2 )~ 3 / 2 , which differs from unity by only 
0.0054 for a near-maximum value of e = 0.06, 
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DIFFERENTIAL PRECESSION: INERTIAL AND DISSIPATIVE COUPLING OF 

THE MANTLE AND CORE 


Statement of Problem 


The hydrostatic figure of a planet represents a compromise between gravitation, which attempts 
to attain spherical symmetry, and rotation, which prefers cylindrical symmetry. Due to its higher 
mean density, the core of the Earth is more nearly spherical than the mantle. The direct lum-solar 
precessional torques on the core will thus be inadequate to make it precess at the same rate as the 
mantle. In fact, the core oblateness is only about 3/4 that required for coprecession with the mantle 
(Smith and Dahlen, 1981). However, it is clearly the case that the core and mantle precess at very 
nearly the same rate (Stacey, 1973). A variety of different physical mechanisms contribute to the 
torques which achieve this coupling, but a purely phenomenological partitioning is useful. The net 
torque can be described as a sum of inertial torques, which are parallel to (Xm x Xc), and dissipative 
torques, which are parallel to (\ m - Xc)- Here, Xc and Xm are the rotation vectors of the core and 
mantle, respectively. The two types of torques have qualitatively different results: inertial torques 
cause the core and mantle axes to precess at fixed angular separations and on opposite side of their 
combined angular momentum vector, whereas the effect of dissipative torques is to reduce the angle 
between the axes. 

On short time scales it is appropriate to consider the core to be an inviscid fluid constrained 
to move within the ellipsoidal region bounded by the rigid mantle (Poincare, 1910; Toomre, 1966; 
Voorhies, 1991). The inertial coupling provided by this mechanism is effective whenever the ellipticity 
of the container exceeds the ratio of the precessional to rotational rates. If the mantle were actually 
rigid, or even elastic (Merriam, 1988; Smylieet al., 1990), this would be an extremely effective type 
of coupling. However, on sufficiently long time scales, the mantle will deform viscously and can 
accommodate the motions of the core fluid (Wu, 1990). The inertial coupling torque exerted by the 
core on the mantle will have the form 

Ti = ki[Xm X Xc] ( 16 ) 

A fundamentally different type of coupling is provided by electromagnetic or viscous torques 
(Rochester, 1962; Sasao et al., 1977; Kubo, 1979). The dissipative coupling torque exerted by the 
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The research design for this segment of the proposed investigation would consist of several parts. 
The objectives and methods are very similar to those described in the previous section. In this 
case, however, the intent would be to determine the precessional and paleoclimatic consequences of 
non-rigd core-mantle coupling, and to explore the possibility that useful constraints on the coupling 
parameters can be obtained from paleoclimatic proxy data. 

A number of estimates already exist for the strengths of inertial and dissipative coupling torques 
(Toomre, 1966; Roberts, 1972; Rochester, 1976; Loper, 1975; Stix, 1982). By most accounts, the 
inertial torque is ~ 5 10 20 N m, and the various viscous and magnetic dissipative torques are 102-104 
times weaker. However, the inertial torque estimates are simply based on the premise that the core 
must coprecess with the mantle. 

Solutions to the coupled precession problem can be found in a form analogous to Equation (11) 
for the rigid precession problem. The principle difference in the present situation is the increased 
richness of the free and forced oscillation spectra. There are modes in which the core and mantle 
precess together, and other modes which reflect differential precession. It is clear that the most 
climatically relevant behavior is the precessional motion of the mantle. Thus the chief interest, from 
that perspective, will be to explore the behavior of the mantle precession modes over plausible range 
of parameter values. Within the geophysically prescribed range of parameter space, are there any 
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climatically significant perturbations to obliquity or longitude of perihelion occasioned by partially 
decoupling the core from the mantle? 

It is very clear that the mantle and core exhibit differential motions on nutational time scales 
(Mathews et ah, 1991; Mathews and Shapiro, 1992). However, one of the difficulties in constructing 
a differential precession model which is truly useful for paleoclimatic studies is that the precise 
geodetic techniques, which are able to constrain nutation amplitudes within a few milliarcseconds, 
still have short enough time spans that most of the differential precession modes of interest are still 
indistinguishable from rigid rotation. 
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DEGLACIAL POLAR MOTION AND ABRUPT CLIMATE CHANGE 


Statement of Problem 

Recent models of the surficial mass transport associated with the last deglaciation (Tushingham 
and Peltier, 1991; Nakada and Lambeck, 1988) suggest that the magnitude of the flow and its 
departure from axial symmetry were both great enough that, if it were not closely balanced by 
internal flow and deformation, it would cause a shift in the body- fixed location of the axis of greatest 
inertia, by an amount of order lo. If the pole were to move by a significant fraction of a degree, the 
resulting displacement of the geographic pattern of land and water relative to the incident radiation 
pattern would cause a climatic perturbation which could either augment or retard the progress of 
deglaciation, depending on the spatio-temporal pattern of the perturbation. The potential thus 
exists for a previously unexplored feed-back loop in the climate system. 

During the peak of the last deglaciation, there was a brief but significant return to full glacial 
conditions. This Younger Dryas climatic event is perhaps best documented in the North Atlantic 
(Ruddiman and McIntyre, 1981; Broecker et ah, 1988), but appears to have been global in scale 
(Currey, 1990; Engstrom et ah, 1990; Gasse et ah, 1991). Though the broad scale timing of the 
deglaciation is consistent with astronomical forcing, the Younger Dryas perturbation was too brief to 
be a linear response to the classical orbital or rotational variations. However, neither the magnitude 
nor duration is a priori inconsistent with a response to deglacial polar motion. The objective of this 
portion of the proposed investigation would be to examine the geodynamic and climatic consequences 
of deglacial mass flow. 

Approach 

The usual geodynamic modeling approach to deglacial polar motion is to use historical observations 
of the rate and direction of polar motion during this century, in conjunction with estimates of the 
surficial mass flow, to derive constraints on deep Earth structure models. However, these same 
models, once they have been calibrated, can be made to deliver estimates of the rates and directions 
of polar motion which occurred during the deglaciation. The magnitude of deglacial polar motion 
depends rather sensitively on the rate and spatial pattern of surficial mass transport and on the 
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internal structure of the Earth, as reflected in the spectrum of relaxation times for surficial loads 
and body forces (Sabadini et al., 1982; Wu and Peltier, 1984; Spada et al., 1992; Ivins et ah, 1992). 

A class of climate models which is very well suited to the proposed investigations incorporates 
spatial and temporal variations in a number of components of the global energy balance. The basic 
equation which governs the energy balance is (North et ah, 1981) 

A + BT + C^-V -{DVT) = QaF (20) 

The first two terms on the left parameterize outgoing radiation, the third term represents storage 
of heat and the fourth term represents divergence of the heat flux. The right hand side is the amount 
of incident radiation that gets absorbed. 

Values for the outgoing radiation parameters (A = 210 W m -2 , B = 2.1 W m -2 K -1 ; Short 
et ah, 1984) can be estimated from satellite data. The heat capacity C is large over water and 
small over land (C^/B = 4.6 years, C / = C w / 60; North and Coakley, 1979). Though much of the 
lateral transport of heat is accomplished by winds and ocean currents, these effects can be modeled 
by a diffusive transport, with D/B = 0.310 (North, 1975; Wyant et ah, 1988). For the incident 
radiation, Q = S/4 = 342 W m -2 (Schatten and Orosz, 1990), F is the spatial and temporal pattern 
of projected area, and a (0.75, Stephens et ah, 1981) is the co-albedo. 

Solutions to this equation with specified forcing can be readily obtained in the spectral domain, 
with spatial patterns represented by spherical harmonic series and temporal patterns represented by 
Fourier series. If A, B and D are assumed to be constants (no spatial or temporal variations), the 
base state solution to Equation (20) is simply 

T nmp = W~ 1 [QG a nmp F n0p ) (21) 

where G (whose 9 indices have been suppressed for typographical clarity) is a coupling constant for 
products of spherical harmonics (Rotenberg et ah, 1959; Dill, 1991) and 

W = [B + n(n + 1) D + ipGC nmp ] (22) 

In the same notation, the temperature perturbation ST due to a small change in radiative forcing 
6F (associated with e, e, or w) is 

6T n 

mp — W~ l [QG dump 6F n 0p ] (23) 

Similarly, the temperature perturbation due to a change <5C in the distribution of land and water 
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relative to the rotation axis (associated with deglacial polar motion) is 

6T nmp = W~ l [—ipG 6C nmp T nmp ] (24) 

Specific questions which should be addressed include: 

• What is the likely history (magnitude and direction) of polar motion induced during deglacia- 
tion? 

• How does the climatic impact of that polar motion compare in magnitude and spatio-temporal 
pattern with the impact of a 1° change in obliquity, or a 1% change in orbital eccentricity? 

• Is the timing and magnitude of the climatic impact such that it would significantly perturb 
the deglaciation? 

• Is this feed-back loop a viable contributor to the Younger Dryas climatic event? 
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FOURIER, LEGENDRE, AND MILANKOVITCH: ESTIMATING DIAGNOSTIC 
SEASONS AND LATITUDES FROM CLIMATE PROXY RECORDS 


Statement of Problem 

As was mentioned in the introduction, a weakness of previous paleoclimate models is the com- 
mon practice of using an overly simplistic representation of the insolation forcing (values for July 
at 65° N are a particular favorite) when attempting to reconcile computed variations in orbital and 
rotational parameters with observed climate proxy records (Hays et al., 1976). An alternative ap- 
proach, which shares many of the same problems is to use a linear combination of the orbital and 
rotational parameters as input to the models (Imbrie and Imbrie, 1980). A problem in using either 
a linear combination of orbital and rotational parameters, or a localized spatio-temporal measure 
of insolation, as a proxy input to climate models is the difficulty in choosing which single value to 
use (Broecker, 1966; Broecker and Van Donk, 1970). Furthermore, the two alternative formats for 
making the choice (time and location versus orbital and rotational parameter combination) are not 
equivalent. It is clear that choice of a specific latitude and time of year to represent the insola- 
tion pattern implies a unique (though non-linear) combination of orbital and rotational parameters. 
However, the converse is not necessarily true. For some combinations of (c,e,zu) there is no corre- 
sponding location and time of year, at which the insolation pattern will be dominated by the selected 
combination of parameters. 

In comparing calculated insolation pattern variations with an observed time series of some climate 
proxy (marine sediment oxygen isotope variations, for example), it is useful to have a simple physical 
model in mind. The observed isotopic anomalies are, to first order, due to two effects; global ice 
volume fluctuations and local to regional scale temperature variations (Emiliani and Shackleton, 
19/4; Kahn et ah, 1981). These effects differ in two important ways. The ice volume effect is 
global in scale and its first time derivative should be proportional to insolation driven temperature 
changes. In contrast, the direct isotopic temperature effect is local to regional in scale, and is directly 
proportional to insolation driven temperature changes. 
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Approach 


A modeling strategy which addresses several of these points involves estimating the linear com- 
bination of Fourier-Legendre coefficients of the insolation pattern which best reproduces the data 
records. This approach provides the opportunity to partition observed variations into global and 
regional effects. It also provides information on which aspects of the insolation pattern variations 
are most diagnostic of the observed proxy variations. If July insolation at 65° N is truly signifi- 
cant- in a particular data record, the selected amplitudes for the Fourier-Legendre coefficients should 
clearly reflect that fact. Alternatively, if the amplitude of the seasonal cycle in the polar regions or 
tropics, or the annual mean equator to pole insolation contrast, or any other linear functional of the 
insolation pattern is most diagnostic, the analysis will indicate that fact. 

As an example of a situation in which this approach could be used, Park and Maasch (1992) 
have recently compared the climate proxy record provided by 6 18 0 records from two long, high 
resolution cores (ODP 677 from the eastern equatorial Pacific (Shackleton et al., 1990) and DSDP 
607 from the mid-latitude North Atlantic (Ruddiman et al., 1989; Raymo et al., 1989) with a 
new estimate of insolation at 65° N (Berger and Loutre, 1991). The three data records analyzed 
(benthic and planktonic records at ODP 677, benthic only at DSDP 607) have much in common, 
but there are also interesting, and possibly significant, differences. It is clear that variations which 
are common to several locations are more likely diagnostic of global climatic variations. Differences 
between locations, or differences between benthic and planktonic records at a single location provide a 
different and complementary view of the climatic response to insolation forcing. Figure (5) illustrates 
a 2 10 6 year history of variations in the amplitudes of the Fourier-Legendre coefficients and Figure 
(6) compares the DSDP 607 and ODP 677 data with a standard insolation curve. 

Several specific questions should be addressed. The first category of questions relate to calibration 
and characterization of technique. 

• Do the Fourier-Legendre coefficient time series form an orthogonal basis? 

• What output time series would correspond to a spatio-temporal delta function (summer solstice 
at 65° N, for example)? 

• What is the resolution (A/i,AM) versus data series length? 

• How does the resolution depend on other factors besides the data string length? 
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• Are two short data strings from significantly different locations better than a single data string 
of equivalent aggregate length? 

The next level of question pertains to the physical meaning or significance of a particular linear 
combination of Fourier-Legendre coefficients. Any such linear combination is equivalent to a spatio- 
temporal pattern of some sort. The result obtained with a finite resolution data structure will 
represent the convolution of the actual physical response with the data resolution kernel. How 
do errors in data chronologies map into errors in resolved spatio-temporal patterns of diagnostic 
insolation? 

The final level of question relates to understanding the implications and significance of results 
obtained by this algorithm from actual data series. 

• Using a variety of oceanic and continental paleoclimatic data series, what are the resolved 
patterns of climatically diagnostic insolation? 

• What physical mechanisms are implicated? 

• Are significantly different patterns obtained from data strings of different age? 

• Is there a single pattern which duplicates the observed increase in amplitude of the 100 kyr 
signal after the Brunhes-Matuyama magnetic reversal? 
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Figure 1 

Orbital eccentricity spectrum and time series. 

Both are based on the secular variation model of Laskar (1988). 
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Figure 2 

Orbital inclination spectrum and time series. 

Both are based on the secular variation model of Laskar (1988). 
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Figure 3 

Rotational obliquity spectrum and time series. 

Time series is based on rotational theory of Kinoshita (1977) and orbital model of Laskar (1988). 
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Figure 4 

Latitudinal and Seasonal Insolation Pattern. 
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Figure 5a 

Insolation Pattern Fourier-Legendre Coefficient Time Series. 
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Figure 5b 

Insolation Pattern Fourier-Legendre Coefficient Time Series. 
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Figure 6 

Paleoclimate Proxy Records and Insolation Time Series. 

Curves are (top to bottom): DSDP 607 (benthic), ODP 677 (planktonic), ODP 677 (benthic), 65° N insolation 
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